
import numpy as np
import pandas as pd
from math import pi,exp

import matplotlib.pyplot as plt
%matplotlib inline


trainset = pd.read_csv('scores_train_2.csv',names=['x1','x2','y'])
testset = pd.read_csv('scores_test_2.csv',names=['x1','x2','y'])
testset.head()


x1 = np.array(trainset['x1'].values)
x2 = np.array(trainset['x2'].values)
m = x1.size
x_train0 = np.concatenate([x1,x2],axis=0).reshape(2,m)
y_train = np.array(trainset['y'].values).reshape(1,m)

x1 = np.array(testset['x1'].values)
x2 = np.array(testset['x2'].values)
x_test0 = np.concatenate([x1,x2],axis=0).reshape(2,m)
y_test = np.array(testset['y'].values).reshape(1,m)

print("Training Set: ", x_train0.shape, y_train.shape)
print("Test Set:     ", x_test0.shape, y_test.shape)

Training Set:  (2, 100) (1, 100)
Test Set:      (2, 100) (1, 100)


def plot_dataset(x, y):
    n = x.shape[1]
    x1 = x[0,:].reshape(1,n)
    x2 = x[1,:].reshape(1,n)
    plot_data(x1,x2,y)
    
    
def plot_data(x1,x2,y):
    indices_pass = np.where(y==1)
    indices_fail = np.where(y==0)

    x1_pass = x1[indices_pass]
    x2_pass = x2[indices_pass]
    x1_fail = x1[indices_fail]
    x2_fail = x2[indices_fail]
    plt.scatter(x1_fail,x2_fail,marker='o',color='green',label='not admitted')
    plt.scatter(x1_pass,x2_pass,marker='x',color='red',label='admitted')
    plt.xlabel('Exam 1 score $x_1$')
    plt.ylabel('Exam 2 score $x_2$')
    plt.legend(bbox_to_anchor=(1.1, 1))
    axes = plt.gca()
    plt.show()


plot_dataset(x_train0, y_train)
plot_dataset(x_test0, y_test)


def normalize(X):
    ### START YOUR CODE
    mu = X.mean()
    std = X.std()
    
    X_norm = (X-mu)/std
    
    return X_norm, mu, std
    ### END YOUR CODE


x_train,mu,stdev = normalize(x_train0)
x_test = (x_test0-mu)/stdev


def dummy_predictor(x):
    rnd = np.random.uniform(size=(1,x.shape[1]))
    return np.round(rnd)


def pass_rate(x):
    n_pass = np.sum(dummy_predictor(x_train))
    rate = n_pass/x_train.shape[1]
    return rate


print(pass_rate(x_train.shape[1]))

0.57


performance_check = [pass_rate(x_train) for i in range(1000)]
print(np.mean(performance_check))

0.49917


def polynomial_features(x, order):
    """
    Arguments:
    x -- input data as numpy array of shape (2,m) where m is the number of samples
    order -- the max order of terms to be added (x1^j*x2^i and i+j<=order)
    
    Returns:
    numpy array of shape (n,m) where n = (order+1)*(order+2)/2 (all the monomials x1^j*x2^i and i+j<=order)
    """
    m = x.shape[1]
    x1,x2 = x[0,:].reshape(1,m),x[1,:].reshape(1,m)
    features = np.concatenate([np.ones((m),dtype='float').reshape(1,m),x1,x2]).reshape(3,m)
    n = 3
    if order > 1:
        for i in range(2,order+1):
            for term in range(i+1):
                features = np.append(features, (x1**(i-term)*x2**term).reshape(1,m), axis=0)
                n += 1
    return features


def initialize_weights(n):
    return np.random.normal(size=(1,n))*0.01


# Auxiliary methods suited for performing the optimize-step below

def predict(X,W):
    """
    Computes the predicted value - given the inpute feature matrix of shape (n,m) and weights vector of shape (1,n).
    The number of features n also includes the bias term.
    """
    ### START YOUR CODE
    h = W.dot(X)
    out = 1/ (1+np.exp(h))
    
    return out
    ### END YOUR CODE
    
def cost(A,Y):
    """
    Computes cross-entropy cost.
    
    Arguments:
    A -- Activations
    Y -- Labels
    """
    m = Y.shape[1]
    c = -np.sum(Y*np.log(A)+(1-Y)*np.log(1-A))/m
    
    return c

def error_rate(A, Y):
    """
    Computes the error rate.

    Arguments:
    A -- Activations
    Y -- Labels
    """
    Ypred = np.round(A)
    
    return np.sum(Y != Ypred) / Y.size

def gradient_cost(A,X,Y):
    """
    Computes the gradient for the cost with respect to the weights vector of size (1,n)
    
    Arguments:
    A -- Activations of shape (1,m)
    X -- Input of shape (n,m)
    Y -- Labels of shape (1,m)
    
    Returns:
    Vector of shape (n,m)
    """
    ### START YOUR CODE
    grad = (Y-A).dot(X.T)
    
    return grad
    ### END YOUR CODE


def optimize(Xtrain,Ytrain,Xtest,Ytest,nepochs,learningrate):
    """
    Implements (batch) gradient descent for minimizing cross-entropy cost. It returns the learning curves 
    for cost and error rate (test and training). The curves are returned as numpy array of lenghth nepochs+1 
    (the +1 for the initial values).  
    
    Arguments:
    Xtrain -- input data for training, numpy array with shape (n,m)
    Ytrain -- labels for training, numpy array with shape (1,m)
    Xtest -- input data for test, numpy array with shape (n,m)
    Ytest -- labels for test, numpy array with shape (1,m)
    nepochs -- number of epochs
    learningrate -- learning rate
    
    Returns:
    traincosts -- learning curve with the cost on the training dataset, a numpy array of shape (nepochs+1)
    testcosts -- learning curve with the cost on the test dataset, a numpy array of shape (nepochs+1) 
    trainerror -- learning curve with the error rate on the training dataset, a numpy array of shape (nepochs+1)
    testerror -- learning curve with the error rate on the test dataset, a numpy array of shape (nepochs+1)
    W -- parameter vector, a numpy array of shape (1,n+1)
    """
    ### START YOUR CODE
    n = Xtrain.shape[0]    
    W = initialize_weights(n)
    
    train_costs = []
    train_errors = []
    test_costs = []
    test_errors = []
    
    for i in range(nepochs):
        # predict the train set
        pred_train = predict(Xtrain, W)
        # wheight update
        W -= learningrate * gradient_cost(pred_train, Xtrain, Ytrain)
        
        # train
        train_costs.append(cost(pred_train, Ytrain))
        train_errors.append(error_rate(pred_train, Ytrain))
        
        # test
        pred_test = predict(Xtest, W)
        test_costs.append(cost(pred_test, Ytest))
        test_errors.append(error_rate(pred_test, Ytest))
        
    return np.array(train_costs), np.array(test_costs), np.array(train_errors), np.array(test_errors), W
    
    ### END YOUR CODE


def evaluate(x_train,y_train,x_test,y_test,pol_degree, nepochs, learningrate):
    """
    Evaluate a model by training it, plotting the learning curves and the decision boundary and 
    returning the performance (final cost and error rate obtained for training and test set)
    """
    Xtrain = polynomial_features(x_train, pol_degree)
    Xtest = polynomial_features(x_test,pol_degree)
    traincosts, testcosts, trainerror, testerror, W = optimize(Xtrain,y_train,Xtest,y_test,nepochs,learningrate)
    plot_curves(traincosts, testcosts, trainerror, testerror)
    Jtrain, Jtest, etrain, etest = traincosts[-1],testcosts[-1],trainerror[-1],testerror[-1]
    print(Jtrain, Jtest, etrain, etest)
    print(W)
    return Jtrain, Jtest, etrain, etest, W


def plot_curves(traincosts, testcosts, trainerror, testerror):
    iterations = range(traincosts.size)
    f = plt.figure(figsize=(10,3))
    plt.subplot(1,2,1)
    plt.plot(iterations, traincosts,label="train")
    plt.plot(iterations, testcosts, label="test")
    plt.xlabel('Epochs')
    plt.ylabel('Cost')
    plt.legend()
    plt.subplot(1,2,2)
    plt.plot(iterations, trainerror, label="train")
    plt.plot(iterations, testerror, label="test")
    plt.ylabel('Error')
    plt.xlabel('Epochs')
    plt.show()


def plot_decision_boundary(x, y, W, pol_degree):
    x1 = x[0,:]
    x2 = x[1,:]
    indices_pass = np.where(y[0,:]==1)
    indices_fail = np.where(y[0,:]==0)

    x1_pass = x1[indices_pass]
    x2_pass = x2[indices_pass]
    x1_fail = x1[indices_fail]
    x2_fail = x2[indices_fail]

    f, ax = plt.subplots(figsize=(7, 7))
    ax.scatter(x1_fail,x2_fail,marker='o',color='green',label='not admitted')
    ax.scatter(x1_pass,x2_pass,marker='x',color='red',label='admitted')
    plt.xlabel('Exam 1 score $x_1$')
    plt.ylabel('Exam 2 score $x_2$')
    ax.legend(bbox_to_anchor=(1.1, 1))
    axes = plt.gca()

    x1_min, x1_max = x1.min() - 1, x1.max() + 1
    x2_min, x2_max = x2.min() - 1, x2.max() + 1
    xx1, xx2 = np.meshgrid(np.arange(x1_min, x1_max, (x1_max-x1_min)/100), np.arange(x2_min, x2_max, (x2_max-x2_min)/100))
    xx10 = xx1.reshape(1,xx1.size)
    xx20 = xx2.reshape(1,xx2.size)
    xx = np.concatenate((xx10, xx20), axis=0) 
    yy = W.dot(polynomial_features(xx, pol_degree)).reshape(xx1.shape)
    
    ax.contour(xx1, xx2, yy, levels=[0], cmap=plt.cm.Paired)
    #ax.axis('off')

    plt.show()


import itertools

pol_degree = list(range(1, 20, 2))
nepochs = [500, 1000]
learningrate = [0.01, 0.1, 0.2]

for pol_d, n_epochs, lr in itertools.product(pol_degree, nepochs, learningrate):
    traincost1, testcost1, trainerror1, testerror1, W1 = evaluate(x_train, y_train, x_test, y_test,
                                                                  pol_d, n_epochs, lr)
    plot_decision_boundary(x_train, y_train, W1, pol_d)
    plot_decision_boundary(x_test, y_test, W1, pol_d)

0.2035284039304025 0.2150584103716009 0.11 0.11
[[-1.68739413 -3.83697297 -3.74561364]]

0.20349770158944003 0.2148453653539708 0.11 0.11
[[-1.7198245  -3.90445628 -3.81433603]]

0.20349770158944003 0.21484536535397059 0.11 0.11
[[-1.7198245  -3.90445628 -3.81433603]]

0.20349778911011124 0.21485515996671223 0.11 0.11
[[-1.71808132 -3.90082616 -3.8106392 ]]

0.20349770158944008 0.2148453653539705 0.11 0.11
[[-1.7198245  -3.90445628 -3.81433603]]

0.20349770158944003 0.21484536535397059 0.11 0.11
[[-1.7198245  -3.90445628 -3.81433603]]

0.07379648134763711 0.13235408493916367 0.0 0.08
[[-5.10391065 -3.8739093  -3.16273569  2.42927413  1.49504784  1.90395244
  -1.38188536 -1.75577284 -1.79873164 -1.96903358]]

<ipython-input-14-2bfc2a07c69d>:24: RuntimeWarning: divide by zero encountered in log
  c = -np.sum(Y*np.log(A)+(1-Y)*np.log(1-A))/m
<ipython-input-14-2bfc2a07c69d>:24: RuntimeWarning: invalid value encountered in multiply
  c = -np.sum(Y*np.log(A)+(1-Y)*np.log(1-A))/m

nan nan 0.0 0.06
[[-15.56262343 -12.33554989 -11.15873562   7.43953941   5.09775387
    6.00607155  -3.18231357  -6.23982521  -6.52618125  -4.80099345]]

<ipython-input-14-2bfc2a07c69d>:24: RuntimeWarning: divide by zero encountered in log
  c = -np.sum(Y*np.log(A)+(1-Y)*np.log(1-A))/m
<ipython-input-14-2bfc2a07c69d>:24: RuntimeWarning: invalid value encountered in multiply
  c = -np.sum(Y*np.log(A)+(1-Y)*np.log(1-A))/m

	x1	x2	y
0	39.196334	78.530294	0
1	40.448499	86.839470	1
2	65.571920	44.303497	0
3	79.648113	70.806564	1
4	66.260221	41.672703	0

Optimizing Model Complexity¶

Load the Data¶

Normalize the Data¶

Dummy Predictor¶

Modelling¶

Predict, Optimize¶

Evaluate the different models¶